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Abstract 

An efficient conditioning technique, the so-called Brownian Bridge simulation, has pre- 
viously been applied to eliminate pricing bias that arises in applications of the standard 
discrete-time Monte Carlo method to evaluate options written on the continuous-time 
extrema of an underlying asset. It is based on the simple and easy to implement ana- 
lytic formulas for the distribution of one-dimensional Brownian Bridge extremes. This 
paper extends the technique to the valuation of multi-asset options with knock-out barri- 
ers imposed for all or some of the underlying assets. We derive formula for the unbiased 
option price estimator based on the joint distribution of the multi-dimensional Brownian 
Bridge dependent extrema. As analytic formulas are not available for the joint distri- 
bution in general, we develop upper and lower biased option price estimators based on 
the distribution of independent extrema and the Frechet lower and upper bounds for the 
unknown distribution. All estimators are simple and easy to implement. They can always 
be used to bind the true value by a confidence interval. Numerical tests indicate that our 
biased estimators converge rapidly to the true option value as the number of time steps 
for the asset path simulation increases in comparison to the estimator based on the stan- 
dard discrete-time method. The convergence rate depends on the correlation and barrier 
structures of the underlying assets. 

Key Words: Monte Carlo simulation, extreme values, Brownian Bridge, multi-asset barrier 
option, multi-variate joint distribution, Frechet bounds. 
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1 Introduction 



Barrier options introduced by Merton (1973) are used widely in trading now. The option 
is extinguished (knocked-out) or activated (knocked-in) when an underlying asset reaches a 
specified level (barrier). A lot of related more complex instruments for example bivariate 
barrier, ladder, step-up or step-down barrier options have become very popular in over-the- 
counter markets. In general, these options can be considered as options with payoff depending 
upon the path extrema of the underlying assets. A variety of closed form solutions for such 
instruments on a single underlying asset have been obtained in the classical Black-Scholcs 
settings of constant volatility, interest rate and barrier level. Sec for example Hcyncn and Kat 
(1994a), Kunitomo and Ikcda (1992), Rubinstein and Reiner (1991). If the barrier option is 
based on two assets then a practical analytical solution can be obtained for some special cases 
considered in Heynen and Kat (1994b) and He, Keirstead and Rebholz (1998). 

In practice, however, numerical methods are used to price the barrier options for a number 
of reasons, for example, if the assumptions of constant volatility and drift are relaxed or payoff 
is too complicated. Numerical schemes such as binomial and trinomial lattices (Hull and White 
(1993), Kat and Verdonk (1995)) or finite difference schemes (Dewynne and Wilmott (1993)) 
can be applied to the problem. However, the implementation of these methods can be difficult. 
Also, if more than two underlying assets are involved in the pricing equation then these methods 
are not practical. 

In this paper we focus on a Monte Carlo simulation method which is a good general pricing 
tool for such instruments. However, finding the extrema of the continuously monitored assets 
by sampling assets at discrete dates, the standard discrete-time Monte Carlo approach, is com- 
putationally expensive as a large number of sampling dates and simulations are required. Loss 
of information about all parts of the continuous-time path between sampling dates introduces 
a substantial bias for the option price. The bias decreases very slowly as 1/\/M for M >> 1, 
where M is the number of equally spaced sampling dates (see Broadie, Glasserman and Kou 
(1997)). Also, extrapolation of the Monte Carlo estimates to the continuous limit is usually 
difficult due to finite sampling errors. For the case of a single underlying asset, it was shown by 
Andersen and Brotherton-Ratchffe (1996) and Beaglehole, Dybvig and Zhou (1997) that the 
bias can be eliminated by a simple conditioning technique, the so-called Brownian Bridge simu- 
lation. The method is based on the simulation of a one-dimensional Brownian bridge extremum 
between the sampled dates according to a simple analytical formula for the distribution of the 
extremum. The technique is very efficient because only one time step is required to simulate the 
asset path and its extremum if the barrier, drift and volatility are constant over the time region. 
We extend the technique to the valuation of multi-asset options with continuously monitored 
knock-out barriers imposed for some or all underlying assets. We derive the general formula 
for the unbiased estimator based on the joint distribution of the multi-dimensional Brownian 
Bridge dependent extrema. In general, however, the analytic formulas are not available for this 
joint distribution and we develop three biased estimators. The upper and lower estimators are 
based on the Prechet bounds for the unknown multi-variate joint distribution of the extrema. 
The third estimator (which is typically most accurate) is based on the joint distribution of 
the independent extrema. The biased estimators can be used to bind the true option price 
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by a confidence interval. Numerical examples indicate that the biases rapidly decrease as M 
increases in comparison to the bias in the standard discrete-time method and the convergence 
rate depends on the correlation and barrier structures of the underlying assets. Finally we dis- 
cuss the application of our biased estimators for the valuation of knock-in, cash-at-hit rebates 
options, lookback options and credit derivatives. 



2 Unbiased estimator via Brownian Bridge correction 
2.1 Model setup 

— * 

Consider the knock-out option Q written on the underlying assets S{t) = {Si{t), Sd{t)). 
The option payout at maturity t = T is some function V{S(T)) if the underlying assets never 
hit the fixed boundaries h(t) < S{t) < H(t), t G [0, T] and zero otherwise. Here H{t) = 
{Hi{t), Hii{t)) and h(t) = {hi(t), hd(t)) are the upper and lower barriers. Hereafter, we 
use vector notation to simultaneously compare all vector components. For example, A < B is 
used to denote Ai < Bi for all i = 1, ...,d. 

Assume that the underlying assets follow risk-neutral geometric Brownian motion 

dSi{t)/Si{t) = fXi{t)dt + ai{t)dWi{t), E[dWi{t)dWj{t)] = Pij{t)dt, (1) 

where S'j(O) is the i-th asset price today, Wi{t),i = 1, ...,d is the ci-dimensional Wiener process, 
/2(t) = (fiiit), fidit)) and (r{t) = {ai(t), ...,adit)) are the drifts and volatilities respectively. 
Let us consider the time slices tm, m = 0, M ordered as = to < < ^2 < ••• < tM = T, 
Stm = tm+i — tm and denote Si{tm) = Sf^\ We assume that the drifts, volatilities, correlation 
coefficients and boundaries are piecewise constant functions of time such that = 
a,{t) = at\ N(t) = Hi{t) = Ht\ Kit) = h^^ if t G [tm, Wi), m = 0, 1, M - 1. Let 
us also introduce the indicator function of the barrier hit at discrete times tm, m = 0, M by 

J _ I 1, if h^""^ < S^^™) < i/^™) for m = 0, M, 
^ \ 0, otherwise. 

In the absence of arbitrage the true option price at t = can be written as an expectation 
Q= f dS^... f rfS(*^)V(S^(*^)MS^ffl|S^(°))-...-p(S("^^ (3) 



where p{S^"^'^^'^\S^'^'^) is the risk-neutral probability density function of the asset value 5^"*+^) 
at tm+i given the asset value S^"^^ at tm- The function should satisfy the Kolmogorov for- 
ward equation (also known as the Fokker-Planck equation) with the absorbing boundaries 
h{t), H{t), t G [0,T], see Cox and Miller (1965). Also we have absorbed the present value 
discount factor into the payoff function V{S{T)) and used the short vector notation for the 
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multi-dimensional Integra^. 

The explicit solution of the Kolmogorov forward equation for the transition probability func- 
tion po('5("+^)|5('")) to unrestricted process ([1]) over the interval [tm, tm+i] (without absorbing 
boundaries) is the (i-variate lognormal distribution with the means: [j^^l^^ — 0.5(crf™'')^](5t„ 



hiSl^^\ variances: {al^^)'^6tm, and linear correlation coefficients: i = l,...,d; j = l,...,d. 

The simulation of 5*1™^^^ from this distribution is simply 



St^ exp [[^^t^ - 0.5{at'>m^ - at^ ^zt^] , (4) 



where Z^^\ i = 1 d are random variates from the (i-variate Normal distribution with the 



linear correlation coefficients pij^\ zero means and unit variances for given m (the random 
variates are independent for different m). 



2.2 Monte Carlo estimators 

The standard discrete-time Monte Carlo approach to estimate the knock-out barrier options ([3]) 
is to assume unrestricted process between sampling dates tm, m = 0, 1, .., M, then simulate N 
independently drawn asset paths according to the iterative equation (j4]) and finally to calculate 
the option price estimate as 

Qs = V{S{T)%^T. (5) 

For simplicity, hereafter, we omit the averaging over N paths. Finding the option price accord- 
ing to (l5l) will introduce a bia^ which is usually larger than the statistical error of the Monte 
Carlo estimates. This is because we lose information on the continuous-time path between 
sampling dates. The bias decreases very slowly as Qs — Q^l/ VM for M >> 1 (see Andersen 
(1996) and Broadie, Glasserman and Kou (1997)). Thus a large number of sampling dates is 
usually required to obtain an accurate estimate of the option price with continuously monitored 
barriers. For example, the bias is still larger than 1% of the true price even for 1024 time steps 
for the case of the standard down-and-out call, see Table 1. Extrapolation to the continuous 
limit is complicated by finite sampling errors of the Monte Carlo estimates. This makes the 
estimator ([5]) computationally expensive. Being mainly interested in eliminating (reducing) the 
biases unaffected by the number of paths, hereafter in formulae, we omit the dependence of 
the Monte Carlo estimates upon assuming that N is large enough to make the statistical 
errors negligibly small in comparison to the biases (we will present the standard errors in the 
numerical examples). 

Consider the continuous-time asset maxima, 

MM = (Mj'"\...,Mi™)), and mmima, 
over the time interval [tm,tm+i], where 

ffj'"' 

1 / = / ds["'\.. J 

/i*"' /i^j'"' 

^ The estimator ([5]) is an unbiased estimator of the option with discretely monitored barriers if sampHng 
dates match the barrier monitoring dates. 
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= max{Siit) : t E [tm, Wi]}, = min{^,(t) : t G [t^, Wi]}- (6) 

If the continuous barriers are imposed between sampling dates then the correct transition 
probabihty function (which is a solution of the Kolmogorov forward equation with absorbing 
boundaries) 

should be used in the option price integral ([3]) rather than the transition probability function 
Po{S^"^~^^^S^'^'>) for an unrestricted process. Once the asset path is simulatecjf] according to (jl]) 
then the unbiased estimator for the knock-out option can be calculated as 

M 

Q = \/(s(T))- w n (8) 

m=l 

where 



is the probability that the assets will not hit the barriers in the time region [tm,tm+i] condi- 
tional on S^"^\ ^C'^+i). This probability is a joint distribution of the extrema L^™') and M^*") 
conditional on S^'^\ S^'^^^\ After suitable normalizatioii^ p(™) is a joint distribution of the 
Brownian Bridge extrema (see Karatzas and Shreve (1991) or Borodin and Salminen (1996)). 
If this probability is easy to calculate then the estimator (j8]) is trivial to implement. One has to 
simulate the asset path at tm, m = 0, 1, M according to the standard procedure (jl]). If the 
underlying assets never hit the barriers at the sampling dates then the option price estimator is 

A/ 

given by the discounted payoff V{S(T)) weig hted with n P^""' and zero otherwise. However, 

m=l 

joint distribution p(™) can be found analytically for some special cases only. In the case of 
barriers discretely monitored at sampling dates, P^"^) = 1 and estimator ([8]) is equivalent to 
the discrete-time estimator ([5]). 



2.3 Marginal distributions of the Brownian Bridge extrema 

For the case of a single barrier for one of the assets at each time region the probability P^*") 
in can be found analytically. Using results from the theory of the Wiener process with 
absorbing boundaries, Cox and Miller (1965), and ([9]) or the formula for the one-dimensional 
Brownian Bridge extremum, Karatzas and Shreve (1991), the marginal distributions of the 

The asset path should not necessarily be simulated according to (|3]). In general it can be simulated from 
any distribution. 

The log-increments of S{t) conditional on 5*^™-' and S'C'^+i) become a so-called Brownian Bridge on 
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maximum M^™^ and minimum of the fc-th asset conditional on S^'^'^ and S''-™^^-' are given 
by 



where 



(10) 
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< min[^f +^4-)], > r^.n[S^r'\st\ 

it\x) = exp (-2 In In -(^/(af))^^^^) . 

\ k k / 

Here ^k^\h^f^^) and ^k^\Hjf!^^) are the marginal probabilities of the upper and lower barrier 
hits by the A;-th asset in the interval [tm,tm+i] respectively. The maximum and minimum 
can be simulated marginally by M^™'-' = — f/) and L^™-* = {^k^^)~'^{U), where ?7 is a 

random variable from the standard Uniform distributioiu. Thus, if the only barrier at [tm, tm+i] 
is imposed for the k-th asset then P^*^) = 1 — ^k^\h^i^^) for the case of a lower barrier and 
p(m) _ for the case of an upper barrier. The marginal distributions (|T0|) are 

valid as long as the interest rate, asset volatility and drift are constant in the time interval 
where the constant barrier is imposed on the asset. 



2.4 Single barrier option 

For a single underlying asset and single barrier per time region it was demonstrated by Ander- 
sen and Brotherton-Ratcliffe (1996) and Beaglehole, Dybvig and Zhou (1997) that simulation 
of the barrier hits in the interval [tm,tm+i] using (ITIl eliminates the bias presented in 
Alternatively we calculate the option price estimator using (jS]). We would like to stress that 
the marginal probabilities ffTOj) can be used to get the unbiased option price estimator not only 
for single asset barrier option but also for multi-asset options if there is a single barrier at each 
time region (this barrier can be imposed for different assets at different time regions). 

In Table 1 we present the Monte Carlo results for down-and-out call for the cases of one 
and two underlying assets. In the first case the option pays max[5'i(T) — K, 0] if Si{t) > 
hi, t G [0,T] and in the second case the option pays max[5'i(T) — K, 0] if S2(t) > /i2, t E 
[0,T]. We have calculated the standard discrete-time biased estimator, Qs, using ([5]) and 
the unbiased estimator, Q, using ([8]) versus the number of equally spaced time steps M. All 
parameters: volatilities, drifts, barriers, correlations are assumed constant. Explicit formulae 
for the unbiased estimators in these examples are 

M 

Q = e-''^max[Si(T) -K,0]- W U ' 

m=l 

for the case of a single underlying asset and 

^ The inverse function (Ci™'')~^(C^) has two solutions. One solution is used to find the maximum and the 
other is used to find the minimum. 
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M 

Q = e-^^max[^i(T) -K,0]- Y[ [1 - (t\h2)] 

m=l 

for the case of two underlying assets. Being mainly interested in eliminating the biases un- 
affected by the number of simulations we did not use any variance reduction technique that 
can be applied to reduce the statistical error of the estimates, see e.g. Boyle, Broadie and 
Glassermann (1997). 

The comparison of the Monte Carlo estimates with the exact option prices calculated by 
analytical formulae demonstrates that Q is an unbiased estimator (the exact value is inside the 
0.95 confidence interval of the estimates for any M). The discretely monitored barrier option 
estimate Qs converges to the continuous barrier case as M increases. However, the convergence 
is very slow and the bias is larger than 1% of the true price even for 1024 time steps. The use 
of the unbiased estimator ([H]) in the above examples is very efficient because one time step is 
enough to obtain the unbiased option price estimate while the standard discrete-time approach 
([5]) requires an enormous number of time steps. 

3 Biased estimators for mult i- barrier option via Frechet 
bounds 

3.1 Frechet bounds for distribution of Brownian Bridge extrema 

The joint distribution P^™) in can be found in closed form via infinite series for the case 
of two dependent extrema (that is two barriers in the same time region) using the results 
from Andersen (1998) for double barrier on a single asset and He, Keirstead and Rebholz 
(1998) for two barriers imposed on different assets. Then, in principle, the unbiased option 
price estimator can be calculated using (jSj). In general, the joint distribution, P*^'"), of three 
and more dependent extrema (that is three and more barriers in the same time region) is 
unknown and unbiased option price estimator ([8]) can not be calculated. However, the univariate 
marginal distributions of the extremes, given by ffTOj) . are known and very simple. The classes 
of multivariate distributions with given margins are the so-called Frechet classes. The results 
on bounds of the distributions with known margins and unknown dependence structure can 
be found in multivariate distribution theory, see for example Joe (1997). The upper and lower 
bounds, so-called Frechet bounds, for the unknown joint distribution with known univariate 
margins are based on simple inequalities involving probabilities of sets. 

Theorem (Lemma 3.8 in Joe (1997)) Let Ai,...,Aj: be the events such that Pr(Aj) = 
a,, i = 1, k. Then 

k 

max[0, ^ai- (k- 1)] < Pr(Ai n ... n A^) < min Oi. (12) 

i=l * 

Let Bi = > /i-™''} and Ai = {Mj^™^ < H^"^^} be the events that minimum is above 

the lower barrier and maximum is below the upper barrier for the i-th asset on [tm,tm+i] 
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respectively. Then Pr(_Bj) = [1 ~ ii^\h!(^'^y\, Pr(y4j) = [1 — ^'f^\H'f^^)\ and the above theorem 
gives the following bounds for the joint probability p(™) in ([9]) 

pM > pM = ^^^[1 _ + it-\h^r\ 0], ^ ^ 

«=i (13) 

pirn) < pM _ _ 1 _ 

j=l,...,d 

The upper bound P^™'' corresponds to a perfect positive dependence between all events Aj, Pj, 
i = 1, If there are only two events then the lower bound P'ff^ corresponds to a perfect 

negative dependence between the events. Perfect positive (negative) dependence between Ai 
and Aj^ Bi and Bj means a perfect positive (negative) dependence between m/™-* and M^'^\ 
llf^^ and L^™-* respectively. Perfect positive (negative) dependence between A^ and Bj means 
a perfect negative (positive) dependence between m/™'' and L^™|^. Also, consider the joint 
probability of the independent event^ Aj, B^, i = 1, d 

= n [1 - ^t^{Ht^)] ■ [1 - ^t^if^)], (14) 

i=l 

which, of course, satisfy Pj^^ < Pj"^^ < pjj^K If all events are positively (negatively) dependent 
then P}™) < P(™) (P)'") > P(")). 



3.2 Frechet bounds and method of images 

The Frechet bounds (fT3l) for the joint distribution P^*") in ([9]) have the following simple inter- 
pretation via the method of images. The joint distribution can be obtained from a solution of 
the Kolmogorov forward equation with absorbing boundaries using The solution for the 
case of a single barrier imposed on the fc-th asset in [tm,tm+i] can be found by the method 
of images, see for example Cox and Miller (1965). The method is based on finding the linear 
combination of the unrestricted process solutions po{S^"'+^^\S^"'^) started at S^^'l the so-called 
source, and Pq{S^'^~^^^ \X'-'^^) started at X^''^\ the so-called primary imag^, satisfying the initial 
and absorbing boundary conditions. The location of the primary image after log-scale change 
is found by reflection of the source in respect to the boundary. This leads to the formula (fTOl) . 
where the contribution of the source is represented by 1 and the contribution of the primary 
image is represented by —^lr\^k^^) (^r ~^fc™''(-f^l™'')) for the case of lower (or upper) barrier. 
Now it is easy to see that the lower bound P^™-* in f |13p is obtained from the source and all its 
primary images (one image for each barrier). The upper bound Pjj^^ in (1131) is obtained from 
the source and one of the primary images that gives the largest contribution. The method of 
images cannot be used to find the exact transition probability in the general case where few 
barriers are imposed on different assets with arbitrary correlation. However, formally, primary 

^The random variables X and Y with continuous marginal distributions have perfect positive (negative) 
dependence if X = T{Y) where r(.) is a strictly increasing (decreasing) function. 

^ This joint distribution has been used by Andersen (1998) and Beaglehole, Dybvig and Zhou (1997) to 
estimate double barrier and double lookback options on a single asset. 

^ In the case of few barriers the primary image may create further images. 
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images can always be introduced and used to get approximate solutions. The numerical results 
below will demonstrate that these approximations are very effective. 

3.3 Three biased estimators for multi-barrier option 

Using the estimators X = L, I ,U defined in ( IT3!) and (fT^ for the joint probability of the 

extremes we can form three biased estimators for the option price ([H]): 

M 

Qx = V{S{T)) ■ I,>T n Px'\ X = U,I, L. (15) 

m=l 

It is easy to see from ([8]), f|T3|) and f|T4l) that 

Ql<Q<Qu, Ql<Qi<Qu (16) 

and, of course, Qx < Qs because P^"^^ < 1, X = UJ,L. While Qi is typically the most 
accurate, Ql and Qu are more useful because they can always be used to bind the true value 
by the confidence interval 

[Ql - ^l~«/2^, Qu + ^l~a/2^], (17) 

where Zi_a/2 is a quantile of the standard Normal distribution, s(.) is the standard deviation, 
N is the number of simulated paths and a is the significance level. Then the true value lies 
within the interval with at least 1 — a probability. The middle of the interval, 

Qo-^Ql + Qu) (18) 

can be used as a point estimator for the true value and half of the interval with -2i_a/2=l can 
be called the standard error of the estimator. Using the inequahty Pj"^ > P(™) (P)"^ < P(™)) 
for negatively (positively) dependent events it is easy to find, for some simple barrier and 
dependence structures, that Qj gives better upper or lower estimator (however, in general, we 
do not know this a priori). Then the true option price can be estimated by 

Qi^UQl + Qi) or Q2^UQi + Qu) (19) 

respectively. The confidence intervals to bind the true value in these cases are formed by 
analogy with (fTTl) . Here we list some examples. If the barrier structure consists of two barriers 
(upper and lower) imposed on one of the assets in each time region then Qi > Q because the 
events of the upper and lower barrier hits are always negatively dependent (that is maximum 
and minimum of the asset are always positively dependent). If the barrier structure consists 
of two upper or two lower barriers imposed on two positively (negatively) correlated assets in 
each time region then Qi < Q {Qi > Q) because the events of the barrier hits are positively 
(negatively) dependent. Also, Qi < Q if only lower or only upper barriers are imposed on 
positively correlated assets in each time region. 
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All three option price biased estimators Qx, X = U, I , L given by flisp are trivial to imple- 
ment because the payoff weights H -^x ? ^ = U,I,L are based on simple marginal distri- 

m=l 

butions ffTOl) . As M increases, it becomes less and less likely that hits of the different barriers 
will occur within the same time interval and the events of barrier hits (and corresponding 
extremes) become disjointed. That is, the probability distribution p(™-) for each of the time 
regions [t^, tm+i] is one of the univariate marginal distributions ([TU]) in the limit tm+i — 0. 
Disjoining of the extrema means that the option price estimator becomes independent from the 
dependence structure (the so-called copula) of the extrema and depends on their marginal dis- 
tributions onljl^. It implies that Qi, Ql and Qu should converge to the true value. Convergence 
of Qi and Qu to each other can be used as a weak criterion of the extreme disjoining. Also 
note that the standard discrete-time biased estimator, Qs, given by (jSD, always converges to 
the true value and Qs > Qu ^ Q- In the next Section we will demonstrate numerically that 
all three estimators Qx, X = U, I, L are rapidly convergent to the true value when compared 
to Q5. The rate of convergence depends on the barrier and asset correlation structures. In this 
paper we do not pursue the analytical derivation of the convergence rate but find it numerically 
for some basic cases. The numerical examples presented in the following Section demonstrate 
the effectiveness of the estimators (fTS!) to correctly price barrier options. 

4 Performance of the biased estimators 

To demonstrate the performance (rapid convergence) of the biased estimators Qx, X = U, I , L 
given in f|T5|) we calculate these estimators and the standard discrete-time estimator Qs, see 
versus the number of equally spaced sampling dates M (the time step is 6t = T/M) for 
the cases of options with two or more knock-out barriers. In addition, we show the results for 
the point estimator Qo ~ (Ql + Qi/)/2, see ( IT8|) . We assume that all parameters (volatilities, 
interest rate, barriers, correlations) are constant and there are no continuous dividends. Being 
interested in the convergence of the biases unaffected by the number of simulations we did not 
use any variance reduction techniques that can be applied to reduce the statistical error of the 
estimates. In all examples the discounted payoff (which is paid at maturity if the assets never 
hit the barriers) is always determined by the first asset: V{Si{T)) = e"''"^ max [5*1 (T) — K,0]. 
The convergence rate should be irrelevant to the payoff paid at maturity. 

4.1 Double knock-out call on a single asset 

First we consider the double knock-out call on a single asset with the lower, hi, and upper. 
Hi, barriers. The explicit expressions for the probability bounds f lT^ . ffT^ used in the biased 
estimators ffT^ are 

^ The phenomenon of maximum and minimum "decouphng" has been noted by Andersen (1998) while using 
independently drawn maximum and minimum to estimate double barrier and double lookback options on a 
single asset. 
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Pf)=max[l-d'"^(/.i)-d"H//i), 0]. 

Actually, the exact joint distribution P^™) is known. It is represented by infinite series (usu- 
ally the series are rapidly convergent) and can be used to calculate the unbiased option price 
estimator, see Andersen (1998). Being mainly interested in the convergence of the biased esti- 
mators we do not pursue this calculation. In Table 2 we show the performance of the tested 
estimators. The maximum and minimum of the asset on the same time interval are always 
positively dependent (the hits of the lower and upper barriers are negatively dependent) thus 
the estimator Qj based on the distribution of the independent extremes is always larger than 
the unbiased estimator Q. In this case Qi is a better upper estimator than Qu. We present the 
results for Qi ^ (Q/ + Ql)/2, see f|T9l) . which can be used as a better point estimator for the 
true price instead of Qq ~ {Qu + Ql)/'^- All our biased estimators Qu, Qi and Ql are rapidly 
convergent to the true value in comparison to the standard estimator Qs- 
Comparison of our biased estimators with the exact analytical result shows that the biases of 
Qi, Ql and Qu are less than their standard errors for M > 4, M > 4 and M > 8 respectively 
while the bias of the standard estimator Qs is significantly larger than its standard error even 
for M=1024. The standard errors are less then 1% of the true value. The exact value is always 
inside the standard confidence intervals of the point estimators Qi and Qq. In the case where a 
double barrier is imposed on the asset the hits of the upper and lower barriers are "physically" 
distant (the asset can not be close to the upper and lower barriers at the same time). Thus 
intuitively we expect that the bias Qu — Ql should be exponentially small for large M. In 
Figure 1 we plot lia{Qu — Ql) versus M for the case of the double knock-out call considered in 
Table 2. The standard error of the plotted estimates is less than the size of the symbols. The 
observed linear behaviour of the graph indicates that the bias decreases as Qu — Ql ^ e~"/^*. 



4.2 Two asset call with two knock-out barriers 

To demonstrate convergence of the estimators for the case where barriers are imposed on 
different arbitrary correlated assets we consider two asset down-and-out call with the lower 
barriers hi and h2 imposed on the first and second assets respectively. In this case the explicit 
expressions for the probabihty bounds (fT51) . ([H]) used in the biased estimators (IT^ are 

Pi-)= min[l-ef)(M], 

j=l,..,a 

P^'^'^=f\[l-iT\h^)], (20) 
pf)=max[l-i:d™H/^.),0], 

i=l 

where d = 2. We designed the problem parameters to equate probabilities of the barrier hits 
^["'\hi) = it\h2) if P = 1. In this case P't'^ = max[l - 2e^^(/ii), 0], P^'"^ = 1 - ^"^(/ii) 
and we expect worst convergence because the events of the barrier hits do not become disjoint 
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at M ^ oo. The exact joint distribution p(™) can be found in a closed form using the results 
in He, Keirstead and Rebholz (1998). It is expressed via infinite series (usually the series are 
rapidly convergent) and can be used to calculate the unbiased estimator ([8]). Again, being only 
interested in the convergence of the biased estimators we do not pursue this calculation. 

In Table 3 we show the exact prices and Monte Carlo estimators for various asset correlation 
values. The exact option price for p = —0.5, 0.5 has been found via numerical integration 
of the two-dimensional density function from He, Keirstead and Rebholz (1998). In the cases 
p = — 1, 0, 1 the option value can be expressed via the barrier options on a single asset 
and the exact solutions are represented via the standard cumulative Normal functional. If the 
correlation, p, between the assets is positive (negative) then the minima of the assets over 
the same time interval are always positively (negatively) dependent. Thus the estimator Qj 
based on the distribution of the independent extremes is always larger (less) than the unbiased 
estimator Q if p < (p > 0). We do not present the results for Qi ^ {Qj + Ql) /2 {li p < Q) 
and Q2 ~ {Qi + Qu)l'^ (if P > 0) but they can be used as a better point estimate for the true 
price instead of Qo ~ {Qu + Ql)/'^- If p = then the minima are independent and Qi is an 
unbiased estimator for the true price because p|"*^ is a valid joint distribution. If p = 1 then 
the minima of the assets have a perfect positive dependence and Qu is an unbiased estimator 
because P^*"^ is a valid joint distribution. 

The obtained results show that our biased estimators Qu, Qi and Ql are rapidly convergent 
to the true value when compared to the standard estimator Qs- Comparison of our results with 
the exact results shows that the biases of Qj, Ql and Qu are less then their standard errors 
if M > 16 while the bias of the standard estimator Qs is significantly larger than its standard 
error even for M=1024 (the standard errors are less then 1% of the true value for almost all 
cases). The only case when the convergence of Qi and is not so rapid is p = 1 (in this 
case Qu is the unbiased estimator). The standard confidence interval of the point estimator Qq 
always contains the exact price. 

The convergence rate of our biased estimators strongly depends on correlation between the 
assets. In Figure 2, Figure 3 and Figure 4 we present the graphs indicating the convergence 
rates for p = — 1, 0, 1. The size of the symbols used for the graphs is larger than the standard 
error of the estimates (we have used more simulations for estimates at large M). 

The linear behaviour of the graphs at large M indicates the following. If p = — 1 then 
Qu — Ql ^ Q-P/^^ for M » 1. This is similar to the results for double barrier call in Table 2 
and Figure 1, because the assets have a perfect negative dependence and the option is equivalent 
to a single asset knock-out option with the flat lower and exponentially growing upper barriers . 
If p = 1 then Qu — Ql^ Vsi for M » 1. This square root convergence is the worst observed 
rate. In this case the asset minima have a perfect positive dependence and the estimator Qu 
is the unbiased estimator of the true option price. That is and Qs converge to the true 
price at the same rate from below and above respectively. If p = then Qu — Ql ~ St"^ for 
M >> 1. We have observed that the convergence rate smoothly deteriorates from the best 

For p = — 1 the option is reduced to a double knock-out call on the first asset with the exponentially 
growing upper and flat lower barriers. For p = the option can be represented as a product of a down-and-out 
call on the first asset and a down-and-out digital option on the second asset. For p — 1 the option is reduced 
to a down-and-out call on the first asset. 
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exponential decay ~ exp{—j3/St) at p = —1 to a rapid power decay St'^ at p = and slow square 
root decreasing ~ ^/6i at p = 1 as the correlation is changed between —1 and 1. Note that we 
have designed the parameters to get worst convergence at p = 1. If we change the parameters 
to make the barrier hit probabilities unequal then the extrema should become disjointed as 
M — > cxo even for p = 1. For example, if we set the asset spots 5*1 (0) = 95, 6*2(0) = 105 and do 
not change the other parameters then the convergence rate at p = 1 is rapid exponential decay 
~ exp(— see Figure 5, while for p = and —1 the rates do not change (that is ~ St"^ and 
~ exp(— respectively). 

4.3 Multi-asset call with multiple barriers 

Finally, to show that our biased estimators work well for real multi-dimensional problems we 
consider, see Table 4, down-and-out call on d underlying assets with the lower barriers imposed 
on all assets for the cases d = 3 and d = 10. The exact (analytical) result is not available 
for this problem. Explicit expressions for the probability bounds (fT3|) . (IT^ used in the biased 
estimators (fTSi) are given by (l20l) . As we have chosen positive correlation between all assets 
the asset minima are positively dependent. Thus Qi is always less than the unbiased estimator 
Q and we present the results for Q2 ~ {Qi + Qu)/'^ which is a better point estimate for the 
true price than Qq. As in the previous examples, all our biased estimators Qu, Qi and Ql 
are rapidly convergent to each other. Their biases become less than the statistical errors for 
M > 16 while the bias of the standard estimator Qs is larger than its standard error even for 
M = 1024. 

5 Conclusions and discussion 

In this paper we have developed a conditioning technique, that can be called a Brownian Bridge 
scheme, for Monte Carlo simulation of a general class of multi-asset options with continuously 
monitored knock-out barriers imposed for some or all underlying assets. We have derived the 
general formula ([8]) for an unbiased estimator of the option based on the joint distribution of 
the multi-dimensional Brownian Bridge extrema If the distribution is known, for example, 
one barrier imposed on one of the assets (or two barriers imposed on different or the same 
assets) per time region, the scheme provides a simple unbiased estimator. The barriers, drifts 
and volatilities are required to be piecewise constant functions of time. In the case of more 
than two barriers per time region the distribution is unknown in a closed form for arbitrary 
dependence between the assets and we derived the upper, Qu, and lower, Qi, biased option 
price estimators. The estimators are based on the Frechet lower, and upper, pIj^\ bounds 

f|T3|) for the unknown joint distribution with given univariate margins. We have also used the 
estimator Qi based on the joint distribution of independent extrema, Pj"^\ For some simple 
barrier and dependence structures Qj can provide better upper or lower bounds. While this 
estimator is usually more accurate, Ql and Qu are often more useful because they can always 
be used to bind the true option price by the confidence interval. As the time between the 
sampling dates decreases, the Brownian Bridge extrema become more disjointed and the biased 
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estimators Ql, Qi, Qu converge to each other and the true value. In the hmit, the option 
price estimator becomes independent of the dependence structure (the so-called copula) of 
the extrema and depends on their marginal distributions only. In our numerical examples we 
showed that the bias Qu — Ql is less and convergence is rapid when compared to the bias and 
convergence of the standard estimator Qs {Qs — Q ~ Vsi for M » 1, where 6t = T/M). 
Usually, Qu — Qu is less than the statistical error of the estimates for only a few time steps. 
In practice, the intermediate dates are introduced due to the interest rate or volatility term- 
structures and the insertion of the additional sampling dates to eliminate the bias may not be 
even required. The convergence rate depends on the barrier and correlation structures. We 
have always observed Qu — Ql ^ exp{—a/6t) for the case of lower and upper barriers imposed 
on the same asset. For the case of two lower barriers imposed on two assets the best detected 
convergence rate is exponential decay ~ exp(— and the worst detected rate is slow square 
root decreasing ~ ^/6i. The worst case was obtained for the unrealistic special set of parameters 
{p = 1 and identical parameters for both assets) which makes barrier hits always equal each 
other and Qu to be an unbiased estimator of the true price. 

The described Brownian Bridge technique is straightforward to use for the valuation of the 
knock-in barrier options and barrier options with constant rebates, R, paid at maturity. The 
unbiased estimators in these cases are given by 

M 

Q = V{S{T))[l-lr>T n ^^"^] 

m=l 

and 

M M 

Q = v(s(T))i.>T n ^^"^ + ^ ■ [1 - ir>T n ^^"^] 

m=l m=l 

respectively. The upper and lower biased estimators can easily be calculated using P^^^ and 
pH) rjj^g technique can easily be applied to efficiently estimate discretely monitored barrier 
options with a large number of observation dates using the method proposed by Andersen 
(1996). That is, using our scheme calculate the continuously monitored barrier option price Qc 
and using the standard Monte Carlo method estimate the option with a low frequency monitored 
barrier. Then the interpolation formula Qm ^ Qc + ^/ VM, where Qm is a barrier option with 
M monitored dates, allows for effective estimation of the option with a high frequency monitored 
barrier. 

Finding the upper and lower biased estimators is not straightforward for the case of multi- 
barrier options with rebates paid at hitting times. This type of problem is also relevant to 
the valuation of credit derivatives. To calculate the unbiased option price estimator, multiple 
hitting times should be simulated from their valid joint distribution which is not known even 
for the case of two barriers. However, the hitting times can easily be simulated from their 
known univariate marginal distributions, see Anderson (1996). Thus, again we have a problem 
of the unknown joint distribution with the known univariate margins. Such options can be 
evaluated in the following way. Calculate one option estimate Q^ simulating the hitting 
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times marginally with perfect positive dependenccj. Another estimate can be found by 
simulating the events marginally with perfect negative dependence in the case of two barriers 
or simulating the events independently if there are more than two barriers. The difference 
between and can be used as a weak criterion for the disjointed events to justify their 
marginal simulations. That is, if the difference is less than the statistical error then the events 
can be assumed disjointed enough (marginal simulation is justified) and {Q^ + Q^)/2 can 
be used as a point estimate for the true price. Otherwise additional sampling dates need to be 
inserted. Lookback type options with payoff dependent on a few continuous extrema can be 
evaluated in a similar way. We will consider these problems in further research. 

In all our examples we have assumed the lognormal diffusion process ([1]) to allow for com- 
parison with analytical results. However, the Brownian Bridge scheme discussed here is still 
applicable for more general diffusion processes. For example, it is applicable if drifts and volatil- 
ities are state-dependent. To solve these models approximate discretisation schemes freezing 
the drifts and volatilities to the left point of each simulation step are used, see for example 
Kloeden and Platen (1992). Then Brownian Bridge schemes can be used because the require- 
ments of piecewise constant drifts and volatilities are satisfied. This will eliminate the bias due 
to discrete underestimation of the continuous extrema (bias due to the discretisation scheme 
will not be removed). 

In this paper we have focused on the application of the technique for the valuation of simple 
knock-out multi-asset options. However, practical use of the technique lies in a broad range of 
problems. For example, the technique can potentially be used for credit and market risk prob- 
lems where valuation of a multi-asset payoff with some barrier levels imposed on the underlying 
assets is very essential. 
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One asset down and out call. 


Two asset down-and-out call. 


The exact value is 8.794. 


The exact value is 8.256 




M 


(5(std.err.) 


(^^(std.err.) 


M 


(5(std.err.) 


(^^(std.err.) 


1 


8.79(0.02) 


10.91(0.02) 


1 


8.26(0.02) 


14.93(0.03) 


2 


8.80(0.02) 


10.66(0.02) 


2 


8.26(0.02) 


13.62(0.03) 


4 


8.80(0.02) 


10.32(0.02) 


4 


8.27(0.02) 


12.35(0.03) 


16 


8.79(0.02) 


9.74(0.02) 


16 


8.27(0.02) 


10.52(0.02) 


64 


8.80(0.02) 


9.33(0.02) 


64 


8.28(0.02) 


9.47(0.02) 


256 


8.80(0.02) 


9.08(0.02) 


256 


8.28(0.02) 


8.90(0.02) 


1024 


8.80(0.02) 


8.94(0.02) 


1024 


8.28(0.02) 


8.59(0.02) 



Table 1: One asset down and out call: 5*1(0) = 100, K = 100, hi = 90, (Ti = 0.3, r = 0.1, 
T = 0.5, 400 000 simulations. Two asset down-and-out call with a single barrier: 6*1(0) = 
^2(0) = 100, K = 100, h2 = 90, r = 0.1, T = 1, ai = a2 = 0.3, p = 0.5, 800 000 simulations. 



M 


Qt/(std.crr.) 


Q/(std.err.) 


QL(std.err.) 


Q5(std.err.) 


(5i(std.err.) 


(5o(std.err.) 


1 


3.01(0.01) 


2.41(0.01) 


1.11(0.01) 


12.23(0.04) 


1.76(0.66) 


2.06(0.96) 


2 


2.21(0.01) 


1.89(0.01) 


1.72(0.01) 


9.60(0.04) 


1.80(0.09) 


1.97(0.26) 


4 


1.84(0.01) 


1.79(0.01) 


1.78(0.01) 


7.41(0.03) 


1.79(0.01) 


1.81(0.04) 


8 


1.79(0.01) 


1.79(0.01) 


1.79(0.01) 


5.73(0.03) 


1.79(0.01) 


1.79(0.01) 


16 


1.78(0.01) 


1.78(0.01) 


1.78(0.01) 


4.50(0.02) 


1.78(0.01) 


1.78(0.01) 


64 


1.78(0.01) 


1.78(0.01) 


1.78(0.01) 


3.06(0.02) 


1.78(0.01) 


1.78(0.01) 


256 


1.78(0.01) 


1.78(0.01) 


1.78(0.01) 


2.40(0.02) 


1.78(0.01) 


1.78(0.01) 


1024 


1.78(0.01) 


1.78(0.01) 


1.78(0.01) 


2.08(0.02) 


1.78(0.01) 


1.78(0.01) 



Table 2: Double knock-out call on a single asset, hi = 900, Hi = 1100, S{0) = X = 1000, 
CTi = 0.2, r = 0.1, T = 0.5, 400 000 simulations. The exact value of the continuously monitored 
option is 1.793. 




-15 I 1 

8 16 ,^ 24 32 



Figure 1: Double knock-out call on a single asset considered in Table 2. 
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p 


M 


(^[/(std.err.) 


(5/(std.err.) 


QL(std.err.) 


(5s(std.err.) 


(5o(std.err.) 


Qc 





1 


5.02(0.03) 


3.65(0.03) 


2.27(0.02) 


11.76(0.07) 


3.64(1.41) 


3.649 




8 


3.78(0.04) 


3.66(0.04) 


3.62(0.04) 


6.84(0.06) 


3.70(0.12) 






16 


3.70(0.04) 


3.66(0.04) 


3.65(0.04) 


5.92(0.05) 


3.67(0.06) 






32 


3.66(0.04) 


3.65(0.04) 


3.65(0.04) 


5.27(0.05) 


3.65(0.05) 






64 


3.65(0.04) 


3.65(0.04) 


3.65(0.04) 


4.81(0.05) 


3.65(0.04) 






1024 


3.64(0.04) 


3.64(0.04) 


3.64(0.04) 


3.93(0.05) 


3.64(0.04) 




0.5 


1 


7.78(0.05) 


5.84(0.04) 


4.22(0.04) 


14.97(0.08) 


6.00(1.82) 


6.527 




8 


6.71(0.06) 


6.48(0.05) 


6.41(0.05) 


10.28(0.07) 


6.56(0.20) 






16 


6.61(0.06) 


6.53(0.06) 


6.51(0.06) 


9.27(0.07) 


6.56(0.10) 






32 


6.57(0.06) 


6.54(0.06) 


6.53(0.06) 


8.52(0.07) 


6.55(0.07) 






64 


6.55(0.06) 


6.54(0.06) 


6.54(0.06) 


7.98(0.06) 


6.55(0.06) 






1024 


6.54(0.06) 


6.54(0.06) 


6.54(0.06) 


6.93(0.06) 


6.54(0.06) 




-0.5 


1 


2.57(0.02) 


1.70(0.01) 


0.67(0.01) 


7.86(0.05) 


1.62(0.96) 


1.395 




8 


1.47(0.02) 


1.41(0.02) 


1.40(0.02) 


3.63(0.04) 


1.43(0.06) 






16 


1.42(0.02) 


1.40(0.02) 


1.40(0.02) 


2.88(0.03) 


1.41(0.03) 






32 


1.41(0.02) 


1.40(0.02) 


1.40(0.02) 


2.45(0.03) 


1.41(0.02) 






64 


1.39(0.02) 


1.39(0.02) 


1.39(0.02) 


2.09(0.03) 


1.39(0.02) 






1024 


1.38(0.02) 


1.38(0.02) 


1.38(0.02) 


1.55(0.03) 


1.38(0.02) 




1.0 


1 


11.36(0.06) 


8.05(0.05) 


6.31(0.05) 


16.79(0.08) 


8.84(2.58) 


11.315 




8 


11.36(0.07) 


10.22(0.07) 


10.00(0.07) 


14.35(0.08) 


10.68(0.74) 






16 


11.37(0.07) 


10.63(0.07) 


10.49(0.07) 


13.63(0.08) 


10.93(0.51) 






32 


11.35(0.07) 


10.84(0.07) 


10.74(0.07) 


13.06(0.07) 


11.04(0.37) 






64 


11.34(0.07) 


10.98(0.07) 


10.91(0.07) 


12.63(0.07) 


11.12(0.29) 






1024 


11.33(0.07) 


11.24(0.07) 


11.22(0.07) 


11.69(0.07) 


11.28(0.12) 




-1.0 


1 


0.415(0.002) 


0.167(0.001) 


0(0) 


2.839(0.018) 


0.207(0.209) 


0.0131 




8 


0.018(0.001) 


0.014(0.001) 


0.014(0.001) 


0.476(0.008) 


0.016(0.003) 






16 


0.014(0.001) 


0.013(0.001) 


0.013(0.001) 


0.250(0.06) 


0.014(0.001) 






32 


0.014(0.001) 


0.014(0.001) 


0.014(0.001) 


0.137(0.004) 


0.014(0.001) 






64 


0.013(0.001) 


0.013(0.001) 


0.013(0.001) 


0.080(0.003) 


0.013(0.001) 






1024 


0.013(0.001) 


0.013(0.001) 


0.013(0.001) 


0.023(0.002) 


0.013(0.001) 





Table 3: Two asset down-and-out call with lower barriers for the first and second assets. 
hi ^ h2 ^ 90, Si{0) = ,^2(0) = 100, K = 100, (7i = (72 = 0.3, r = 0.1, T = 1, 100 000 
simulations. Qc is the exact value of the option with continuously monitored barriers. 
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d 


M 


Qu{std.eTT.) 


(^/(std.err.) 


QL(std.err.) 


Qsistd.evT.) 


(52(std.err.) 


(5o(std.err.) 


3 


1 


8.96(0.07) 


6.69(0.06) 


5.13(0.06) 


14.96(0.10) 


7.83(1.20) 


7.04(1.97) 




2 


8.26(0.07) 


7.20(0.07) 


6.76(0.07) 


13.27(0.09) 


7.73(0.60) 


7.51(0.82) 




4 


7.83(0.07) 


7.43(0.07) 


7.31(0.07) 


11.81(0.09) 


7.63(0.27) 


7.57(0.33) 




8 


7.65(0.07) 


7.51(0.07) 


7.47(0.07) 


10.76(0.09) 


7.58(0.14) 


7.56(0.16) 




16 


7.60(0.08) 


7.56(0.08) 


7.54(0.08) 


9.96(0.09) 


7.58(0.10) 


7.57(0.11) 




32 


7.60(0.08) 


7.59(0.08) 


7.58(0.08) 


9.29(0.08) 


7.59(0.08) 


7.59(0.09) 




64 


7.60(0.08) 


7.59(0.08) 


7.59(0.08) 


8.80(0.08) 


7.59(0.08) 


7.59(0.08) 




1024 


7.60(0.08) 


7.60(0.08) 


7.60(0.08) 


7.91(0.08) 


7.60(0.08) 


7.60(0.08) 


10 


1 


4.62(0.05) 


1.19(0.02) 


0.21(0.01) 


10.36(0.09) 


2.90(1.75) 


2.41(2.23) 




2 


3.56(0.05) 


1.97(0.03) 


1.33(0.03) 


7.92(0.08) 


2.77(0.84) 


2.45(1.16) 




4 


2.98(0.05) 


2.39(0.04) 


2.20(0.04) 


6.13(0.07) 


2.68(0.34) 


2.59(0.44) 




8 


2.80(0.05) 


2.60(0.05) 


2.54(0.05) 


5.09(0.07) 


2.70(0.15) 


2.67(0.18) 




16 


2.71(0.05) 


2.64(0.05) 


2.61(0.05) 


4.37(0.06) 


2.67(0.08) 


2.66(0.10) 




32 


2.67(0.05) 


2.65(0.05) 


2.64(0.05) 


3.84(0.06) 


2.66(0.06) 


2.66(0.07) 




64 


2.65(0.05) 


2.64(0.05) 


2.64(0.05) 


3.48(0.06) 


2.65(0.05) 


2.64(0.06) 




1024 


2.65(0.05) 


2.65(0.05) 


2.65(0.05) 


2.86(0.05) 


2.65(0.05) 


2.65(0.05) 



Table 4: Knock-out calls on three and ten assets with lower barriers for each of the asset. 
Si{0) = 100, hi = 80, pij = 0.5 if i 7^ j, at = 0.4 = 1, ...,d), T = 1, K = 100, r - 0.05, 
100 000 simulations. 




32 64 ^ 96 128 

Figure 2: Two asset down-and-out call considered in Table 3 with p — —1. 
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Figure 3: Two asset down-and-out call considered in Table 3 with p — 0. 




Figure 4: Two asset down-and-out call considered in Table 3 with p — 1. 
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32 64 ,^ 96 128 



Figure 5: Two asset down-and-out call considered in Table 3 with p — 1 and 5*1 (0) = 95, 
^2(0) = 105. 
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